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We consider a system of interacting spin-one atoms in a hexagonal lattice under the presence 
of a synthetic gauge field. Quenching the quadratic Zeeman field is shown to lead to a dynamical 
instability of the edge modes. This, in turn, leads to a spin current along the boundary of the system 
which grows exponentially fast in time following the quench. Tuning the magnitude of the quench 
can be used to selectively populate edge modes of different momenta. Implications of the intrinsic 
symmetries of Hamiltonian on the dynamics are discussed. The results hold for atoms with both 
antiferromagnetic and ferromagnetic interactions. 
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Topological Bloch bands and their concomitant pro¬ 
tected edge modes play important physical roles in sev¬ 
eral solid-state materials including quantum Hall systems 
[1] and topological insulators/superconductors [2, 3]. Re¬ 
cent years have experienced remarkable progress in engi¬ 
neering systems which possess topological Bloch bands 
as a result of induced effective gauge fields. Effective 
gauge fields have been realized in systems of ultracold 
atoms through mechanical rotation, optical lattice ‘shak¬ 
ing’, and laser-assisted tunneling (see [4] for a review). 
Recent milestones include the experimental realization 
of the Rice-Mele [5], Hofstadter [6-8], and Haldane [9] 
models in optical lattice systems. The emerging field of 
topological photonic lattices offers a separate and com¬ 
plementary system where topological Bloch bands have 
also recently been realized [10-13]. Most of the exper¬ 
imental systems with ultracold atoms involve tempera¬ 
tures for which the particle dynamics can be accurately 
described by non-interacting theories. However, as exper¬ 
imental techniques are refined and quantum degeneracy 
is reached, interactions and particle statistics will play 
an essential role in the physics of these systems. 

Of central interest in each of these models is the pres¬ 
ence of topologically protected edge modes. The physical 
consequence of edge modes for a degenerate Fermi gas in 
a topological band is clear: if the Fermi energy resides 
within the bulk band gap, the near-equilibrium dynamics 
will be described by degrees of freedom localized along 
the edge of the system. On the other hand, for bosonic 
systems, particles will generally condense into the lowest 
bulk band, leaving the higher-energy edge states largely 
unimportant for the dynamics of the system. 

In this work, we propose a scheme to bring to the 
fore the role of edge states in the dynamics of bosons 
in a prototypical two-dimensional topological band sys¬ 
tem. In particular, we describe how a quantum quench 
of an interacting spinor generalization of the well-known 
Kane-Mele model [14] can result in an exponentially fast 
growth in the population of edge (but not bulk) states. 
This will be exhibited as an exponential growth of the 


spin current along the boundary of the system. Further¬ 
more, the momenta of these edge states can be selected 
by tuning the magnitude of the quench. 

Before proceeding, we briefly comment on previous 
related work. Probing the topology of band systems 
through a quench has received considerable attention in 
the recent theoretical literature ( e.g . [15-22]). One of 
us [23] has described the exponentially fast population of 
an edge mode in perhaps the simplest topological system, 
the Su-Schrieffer-Heeger (SSH) model [24], by preparing 
a bosonic gas in an excited spatial mode. A photonic 
version of the SSH model was in fact realized [25] where 
on-site absorption was shown to lead to fast population of 
edge states [26]. Quadratic fermionic hamiltonians have 
been previously classified in terms of their symmetries 
[27, 28]. Much less work, however, has been carried out 
for bosonic systems. An expression for the Chern num¬ 
ber, accounting for the symplectic transformation prop¬ 
erties of bosonic systems was obtained in [29]. In [30], 
the topology of real Bogoliubov excitations in inversion- 
symmetric lattices was analyzed. However, to our knowl¬ 
edge, a classification of quadratic bosonic systems allow¬ 
ing for dynamical instabilities does not presently exist. 

The Haldane model [31] was realized in an ultracold 
atomic system by the Esslinger group [9]. This sys¬ 
tem involves atoms on a hexagonal optical lattice where 
spin-orbit coupling is induced through shaking. In this 
work, we consider a spin-one version of this system where 
atoms experience spin-orbit coupling proportional to the 
^-component of their spins. In particular, we introduce 
the following generalization of the Kane-Mele model [14]: 

Usikm = -wJ2 *1 + iX Y, s Aj ( 1 ) 

(ij) « ij » 

where <f>i = (4>^i, is a vector composed of 

bosonic annihilation operators at site i for each spin com¬ 
ponent and we denote the 3x3 spin-one matrices as 
S = (S x , S y , S z ). The second term above describes hop¬ 
ping between second neighbors, and Vij = +1(—1) if 
the atom makes a left (right) turn to reach a second- 
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neighbor site [14]. The spin components are decoupled 
in ^sikm- the spin-zero component is described by the 
nearest-neighbor graphene model with hopping w while 
the spin-±l components are described by two Haldane 
models with opposite magnetic fields. We restrict our 
attention to w > V3\X\ so that the single-particle state 
of lowest energy occurs at the centre of the Brillouin zone. 

Next, we include on-site interactions that preserve spin 
rotation invariance [32-34]: 


H 


“ = £ I (*!*<)' + T (*' s4 ‘) 


( 2 ) 


where U and U s describe the magnitude of the den¬ 
sity and spin interactions respectively. Finally, we in¬ 
troduce the standard quadratic Zeeman effect for spinor 
condensates [34]: 7i ext = gJA (^z) 2 &i- While exter¬ 
nal magnetic fields will provide only positive values of 
g, microwave fields can be utilized to access both pos¬ 
itive and negative quadratic Zeeman shifts [35]. Intro¬ 
ducing a chemical potential /i, the full Hamiltonian reads 
H = HsiKM+Hint+^ext-/iZ^i which is invariant 

under time reversal and global spin rotations about the 
z- axis, the importance of which will be addressed below. 

In this work, we consider a quantum quench that 
abruptly changes the quadratic Zeeman energy q from 
an initially large and positive value to a final value g/. 
This form of quenching has been experimentally achieved 
in several experiments in the past decade (see [34] and 
references therein). The initial state is a coherent state 
with all bosons in a spatially uniform spin-zero state: 

|4r in )=e-3 JV »e VS 2:*Ho|0) (3) 


where N p is the total atom number. There are no par¬ 
ticles with s z = ±1 in this state. This initial state is 
the mean-field superfluid ground state of the graphene- 
lattice boson Hubbard model. A variational calculation 
shows that the average number of bosons per site, n, is 
related to the chemical potential by fi = nU — 3 w. 

To investigate the ensuing dynamics after the quench, 
we consider small fluctuations of the Hamiltonian with 
q = qf around the initial state (3). Let fa = 3b — 
(0,^,0) T where fa = (fai, fa,o, fa^-i) T . The Hamil¬ 
tonian can be expanded to quadratic order in fa as 
H = (^ in |?i|T in ) + H B . One finds 


U b = -w'Y^ 0jq + iX p ij4>\ s z4>j + 

(ij) (W» i 

+ ( -^-fa,o4>i,o + j + H.c. (4) 


where M = diag(3re + I/ s n + g/,3re + I/n,3re + I/ s n + g/) 
and qq is the quadratic Zeeman energy after the quench. 
The ^-rotation symmetry of H B ensures that the spin- 
±1 components are decoupled from the spin-0 compo¬ 
nents and hence the (4) can be written as H B = Ho + 


H± i. The Hamiltonian Ho describes the dynamics of the 
spin-zero components and has no spin-orbit coupling. It 
is readily diagonalized by a Bogoliubov transformation. 
The resulting spectrum is stable and exhibits the usual 
linearly-dispersing phonon mode. From now on, we will 
focus on the spin-±l sector described by H±\. 

As we are primarily interested in the dynamics of the 
edge states of this model, we will focus on the strip geom¬ 
etry. Denoting the primitive lattice vectors of graphene 
as ai and a 2 , we consider open (periodic) boundary con¬ 
ditions along the ai (a 2 ) direction. It is instructive to 
rewrite (4) in the eigenbasis of the non-interacting spin- 
one Kane-Mele model (1) in this geometry. The spin-±l 
Hamiltonian becomes 

^±1 = ^ \i e k ^ — l) 

k,v 

T — l T H.c.)]. (5) 

Here, are the single-particle energies of the Haldane 
model in the strip geometry, k = k • a 2 is the momen¬ 
tum along the periodic direction, and a k:U ,m annihilates 
a boson in the eigenbasis of (1). We have introduced 
A = —U s n — 3 w — qf which serves as the tuning param¬ 
eter for our quench. Due to time-reversal symmetry and 
the spatial uniformity of the initial state, H ±i can be 
separated into pairwise couplings between (fc, zq m) and 
(—fc, i/, — m) modes which greatly simplifies the analysis. 

Unlike Ho, H±i cannot in general be brought to di¬ 
agonal form and may exhibit a dynamical instability. 
Therefore, we focus instead on the Heisenberg equa¬ 
tions of motion: idt&k,v,±i(t) = [&k,v,±i(t),H±i] where 
6tk,v,±i(t) = e^ ±lt d/ CjZy? ±ie _z ^ ±lt and we have set h = 1. 
These can be solved to give 

a k ,v,i(t) — Ak,v{t)oLk,v,i T B k ,zy(^)d_£, v , 

, ^ (6) 
^ — k,u, — 1(^) — -^k,u{l^)^k,u, 1 A ^k,u — k,v, — 1 • 

where A k ^(t) = cos (E^t) - i(e^ - A) sin (E^t)/ 
and B k ^(t) = — iU s nsm(E^t)/E^ with the Bogoli¬ 
ubov energies 

4 V) = ^(^-A) 2 -^) 2 . (7) 

For sufficiently low condensate depletion, the Bogoliubov 
Hamiltonian can be used to propagate the initial state (3) 
as |T(t)) = e _z7 ^ Bt |T in ). The solutions (6) can then be 
used to obtain the expectation value of bilinear operators: 

(4/(t) |d|, | T(t)) = 3k,k'3m,m'\B k ^ u | (8) 

(^Wl & k ,u,m^k' ,u' ,m' l^(^)) — &k,-k'&v,v'&m,-m' A k , v B k , v 

for spin m = ±1 components. 

We now arrive at the central result of this work. If we 
tune the quench parameter A to satisfy 

ejy - u s n < A < gM + u s n, 


(9) 
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FIG. 1. (Color Online) The Bogoliubov energy spectrum 
EE^ corresponding to (7) for the interacting spin-one Kane- 
Mele model in the strip geometry. Parameters in (a) are 
U s n — 0.2 w, A = w/2 and qf — U s n which corresponds to 
a shallow quench with stable spectrum. Parameters in (b) 
are the same as (a) except qf — —3 w — 3 U s n for which the 
bulk states are stable while the edge states experience an ex¬ 
ponentially fast population growth. Gray curves indicate hole 
bands while black curves indicate bulk particle bands. Red 
and blue curves indicate edge states propagating in opposite 
directions. For clarity, edge states on only one side of the 
system are plotted (on the opposite side, the roles of the par¬ 
ticle and hole edge bands are reversed). Imaginary parts of 
eigenvalues are given by dashed lines. 


the Bogoliubov energies E^ becomes imaginary signify¬ 
ing a dynamical instability for the (fc, v) mode. Physi¬ 
cally, this provides exponentially fast population of the 
unstable modes. In other words, through (9), the quench¬ 
ing protocol gives a “window” (centered on A and of 
width 2 U s n) of unstable modes in the spectrum e^\ 

To understand the instability criterion, we now discuss 
the single-particle eigenstates for the spin-ztl compo¬ 
nents. As already mentioned, these are the eigenstates of 
the Haldane model with opposite fluxes for the two spin 
components. The bulk states form two bands of states 
separated by a band gap of 6^/3 A centered about zero 
[31]. In addition, the non-interacting spectrum exhibits 
topologically protected edge modes which exist within 
the bulk gap. In the Supplemental Material, we show 


that the dispersion of the edge states is 

^.edge _ _j_ 6 wAsin (k) 

£ k — ^ /- —— v iU ; 

yjw 2 + 16A 2 sin 2 (|) 

where spin-ztl modes propagate in opposite directions. 
If we tune A to sit inside of the band gap, then, for suffi¬ 
ciently small U s n, we can achieve the intriguing situation 
where the bulk states are stable while the edge states ex¬ 
perience exponentially fast population growth. A similar 
scheme for populating edge modes has been previously 
reported in [23]. Moreover, unlike [23], due to the tune- 
ability of A, the current scheme allows one to selectively 
populate states with particular momenta along the edge. 
From (7), we see that the most unstable modes occur 
at momenta for which = A, so, for instance, when 
A = 0 , edge modes with momenta k = tt will be popu¬ 
lated most rapidly. Bulk and edge bands are shown in 
Fig. 1 for two particular quenches. 

The symmetries (inversion, time reversal, and spin ro¬ 
tation) of the Bogoliubov Hamiltonian greatly simplify 
the above analysis. That is, given the single particle ener¬ 
gies of the Haldane model, the above analysis essentially 
reduced to solving a two-mode problem. For more gen¬ 
eral couplings with less symmetry, it is easiest to proceed 
by solving the Bogoliubov-de Gennes (BdG) equations. 

T 3 H k v k ,„ = E^v ktV (11) 

where H k is the BdG Hamiltonian which, for our prob¬ 
lem, can be directly determined from (4). For a sys¬ 
tem having length 7V|ai| along the direction with open 
boundary conditions, H k is a 2A f x 2A f dimensional ma¬ 
trix where J\f = 6N (the factor of 6 accounts for the spin 
and sublattice degrees of freedom) while 73 = cr 3 0 l/v'x.AG 
The BdG Hamiltonian generically possesses a ‘particle- 
hole’ symmetry which requires the eigenvalues to come 
in pairs. For each pair of stable (real) eigenval¬ 

ues, one member will have positive norm defined with 
the t 3 metric, v k u+ r^v k ^+ > 0 , while the other will have 

negative norm v k ,v~ < 0- In analogy with BCS 

superconductors, we refer to bands composed of eigen¬ 
states with positive (negative) norms as ‘particle’ (‘hole’) 
bands, which are indicated in Fig. 1. 

For stable systems, positive (negative) norm states cor¬ 
respond to positive (negative) eigenvalues E^\ This is 
not the case for unstable systems. In [36], the origin of 
a dynamical instability was traced to positive and nega¬ 
tive norm states that become degenerate in the absence of 
pairing (non-particle number conserving) terms. Then, 
pairing terms generally lift such degeneracies and lead to 
complex Bogoliubov energies. This is precisely the mech¬ 
anism leading to the unstable edge modes in our problem. 
On the other hand, as is evident from Fig. 1(b), one can 
have overlap between bulk particle and bulk hole bands 
which do not lead to dynamical instabilities. 
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We interpret these bulk-band degeneracies as being 
protected by symmetries of the problem. Indeed, in 
the Supplemental Material we show that small contri¬ 
butions to the Hamiltonian that break time-reversal, in¬ 
version, or S z symmetry can hybridize the bulk particle 
and hole bands leading to bulk dynamical instabilities 
which would obscure the population growth of the edge 
modes. Such symmetries can be used to construct other 
bosonic models having unstable topological edge modes 
with stable bulk modes. However, for definiteness and 
due to experimental relevance [9], we focus on the inter¬ 
acting S1KM model in this work. 

We now move on to discuss the physical consequences 
related to the quenching protocol. We first consider the 
number of particles excited into the spin-=bl modes as a 
result of the quench. Using (8), we find 

A4 i(i) =Y,(Vmlmhm\y(t)) =2Y,\Bk,v\ 2 . (12) 

i,m=zLl k,v 


For quenches satisfying (9) and chosen to select only edge 
modes for instability, the population growth in Af±% will 
be localized to the edges of the system. Keeping only 
unstable modes and linearizing the edge spectrum about 
k = 7r, for \U s \nt 1 one finds 


A/±i (t) 



16V n 2 [WM r 2\u.m 
w 2 12A V 7 rt 


(13) 


where N 2 is the number of lattice sites along the a 2 direc¬ 
tion. Note that this expression is also valid for negative 
E/ s , e.g. for 87 Rb atoms. 

The quenching protocol is also expected to create a 
spin current along the edge. The continuity equation for 
the local spin moments S z &i gives an expression for 
the spin current operator. At long wavelengths, the spin 
current operator along the edge is found to be [37] 




where &k,m is a 2 TVi-dimensional vector composed of 
annihilation operators for spin-m bosons on sites in 
a unit cell of the strip geometry. H^ is the non¬ 
interacting matrix Bloch Hamiltonian for spin compo¬ 
nent m which can be directly determined from (1). We 
wish to evaluate the expectation value of this operator 
with the state |T(£)). Writing in an eigenbasis of 

the non-interacting Hamiltonian, employing (8) and the 
Feynman-Heilman relation, we find the intuitive relation 

= (A=o) = J- ^dvetflBk'At)] 2 (15) 

k' ,u 

where the two spin components have contributed an 
equal amount. Additionally, the k ^ 0 components of 


(Jk^) vanish. Under the same conditions as were used 
for the evaluation of A/±i, one finds 

I./jUtU \Us\nt 

4V 7 Tt 

Before closing, we comment on the experimental feasi¬ 
bility of the above protocol. It is expected that quantum 
degeneracy in topological optical lattice systems will be 
reached in the near future with advances in experimental 
techniques. This will open new doors for exploring the 
non-equilibrium dynamics of topological band systems. 
In addition to the optical lattice, atoms are further con¬ 
fined in experiments by an overall (typically harmonic) 
trapping potential. Harmonic traps can obscure the edge 
states in the single-particle spectrum [38] of such systems. 
This problem can be surmounted by using box-shaped 
traps which are now available [39, 40]. On the other 
hand, we expect many aspects of our results with open 
boundary conditions to be qualitatively correct for a har¬ 
monic trap provided the initial state is in the Thomas- 
Fermi regime. The reason is that the condensate will 
screen the trap potential, leading to an effective poten¬ 
tial with sharp boundaries in the BdG equations for the 
edge excitations of interest. Indeed, the dynamics of a 
spinor condensate in a harmonic trap following a quench 
in the quadratic Zeeman field can be accurately mod¬ 
eled by an effective spherical-box potential [41]. Due to 
the sharp boundaries of the effective potential, we ex¬ 
pect that well-defined edge states will continue to exist 
in this geometry. Finally, we note that our results will 
hold for either antiferromagnetic ( U s > 0) or ferromag¬ 
netic ( U s < 0) interactions and so are relevant for both 
87 Rb and 23 Na condensates. 

In summary, we have proposed a method whereby 
topological edge modes are populated exponentially fast 
through a quantum quench. Although edge states are 
typically unimportant for bosonic gases near equilibrium, 
we have shown that the non-equilibrium dynamics after a 
quench can be dominated by degrees of freedom localized 
on the boundary of the system. The growth of the edge 
modes will be limited at longer times by interactions not 
captured in the Bogoliubov theory. The long-time decay 
mechanism of the dynamically populated edge modes will 
also be due to these interaction terms. 

This quenching protocol provides a means of collecting 
cold atoms coherently in a quasi-one-dimensional struc¬ 
ture without the need for extra trapping lasers. Though 
this work focused on ultracold atoms, it will also be 
worthwhile to consider the parallels with photonic lat¬ 
tices where pairing terms can be generated by nonlinear 
optical methods. 
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SUPPLEMENTAL MATERIAL 
Derivation of the edge-state dispersion 

Our aim is to analyze the single-particle edge states of 
^sikm from the main text. Note that since this Hamil¬ 
tonian is invariant under spin rotations about the z axis, 
we can consider each spin component separately. We first 
consider the Hamiltonian describing m = 1 spin compo¬ 
nent, which is the Haldane model. The primitive lat¬ 
tice vectors, as shown in Fig. 2 are ai = (yy, ^ and 

3.2 = where we have set the lattice constant 

to unity. We impose periodic boundary conditions along 
the a 2 -direction and take the system to be semi-infinite 
along the ai-direction. The single particle energies Sk are 
found by solving the following Schrodinger equation 

Rktpn,k + V k *l> n+hk + Vl V’n-l.fc = £ k 1p n ,k (17) 
for n > 1 and 

Rk'tpn,k T Vk'*Pn-\-l,k = ^k^Pn^k ( 18 ) 



FIG. 2. A portion of the hexagonal lattice used for the spin- 
one Kane-Mele Hamiltonian. Triangular sublattice sites A 
and B are labelled by closed and open circles respectively. 
The lattice basis vectors are given by ai and a 2 . 

These equations can be directly solved to determine xpi 
(up to a prefactor), Kq, and e^. In particular, for the 
edge dispersion we find 


for nm 1. Here ip n ^ = (4^2 , Plfk ) T * s a two-component 
vector where n labels unit cell and A , B label sublattice 
(see Fig. 2). Rk and V& are 2x2 matrices given by 


/ — 2Asin(fc) — w(l + e lk ) \ 
y — w(l + e lk ) 2Asin(fc) J 


(19) 


and 
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—iA(l - e ~ ik ) 0 

—w iA(l — e~ lk ) 


( 20 ) 


The Schrodinger equation written above and its con¬ 
stituent matrices directly follow from the non-interacting 
Hamiltonian Hsikm- 

We search for a solution of the form 




1 

^1,k ^2 ,k 


[(«1 , k) n ~ («2,*;)")] 


( 21 ) 


where we require that \n 2 ,k\ < 1 so that the 

solution decays into the bulk. Inserting this into 
the Schrodinger equation, we find (Rk + 14 ^1 ,k + 
^V fe t )-0i,fe = £fet/»i,fe and (Rk + V k K 2 ,k + = 

£k^Pi,k which can be rearranged to give the following two 
equations: 


V k 1V k‘ t P hk = Kl,fcK2,fc^l, (22) 

(Rk + («l,fc + K2,k)V k ) Ipl,k = Zk'lplik- (23) 


£k = ~ 


6w A sin(fc) 

^Jw 2 + 16A 2 sin 2 (|) 


(24) 


The edge mode energies for the m = 0 and m = — 1 spin 
components can be found with the replacements A -> 0 
and A — A respectively. 


Symmetry protection of bulk instabilities 

We now show that small perturbations to the Hamil¬ 
tonian Rb (Eq. 4 from the main text) that break time 
reversal symmetry (TRS), inversion symmetry (IS), or 
S z symmetry can lead to the emergence of dynamical in¬ 
stabilities coming from the hybridization of bulk modes 
in the particle-hole picture. Below we give three terms 
that we add to the Hamiltonian Rb which break either 
TRS, IS, or S z symmetry while preserving the other two. 
To break S z symmetry we consider the contribution of 
same spin-component pairing terms 

Vs z = $ ^2 + <Pi,-l<Pi,-l + H.C.^ (25) 

i 

where S is taken to be a real and small in comparison to 
other energy scales of the problem. This term preserves 
TRS and IS. Inversion symmetry can be broken, while 
keeping the TRS and S z symmetry intact, by adding a 
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FIG. 3. Bogoliubov energy spectrum for the interacting spin-one Kane-Mele model with broken time-reversal, inversion, 
and S z symmetries. Parameters are chosen to be U s n = 0.2 w, A = w/ 2, qf = —3 w — 3U s n/2 , and 5 = 0.095ie. Gray curves 
indicate hole bands while black curves indicate bulk bands. The energy values corresponding to dynamically unstable modes 
(with finite imaginary component) are colored in green. For clarity, modes from the zero spin component (which is always 
stable) are excluded from these plots. 


staggered pairing between the two triangular sublattices 
A and B (Fig. 2): 




(A)UA) 
n . 1 ( Pn ,-1 


Cift+H.c.) (26) 


where n labels unit cells. Finally, we consider a term 
that changes sign under time reversal, but preserves IS 
and S z symmetry: 

v r =i5 ^4>\s 2 z 4>j (27) 

<fe')) 


In Fig. 3 we show the real and imaginary parts of Bo¬ 
goliubov energy spectrum computed numerically after 
adding each of these terms to the Hamiltonian 1~Lb • In 
addition to edge dynamical instabilities discussed in the 
main text, some bulk bands acquire a small imaginary 
part in the energy spectrum. This effect is amplified 
with the increase of 6 and disappears in the absence of 
the symmetry breaking terms. 

We now move on to discuss a slightly more stringent 
sufficient condition for the absence of bulk instabilities. 
We impose the condition of TRS and S z symmetry, and 
further require that the mean-field density (used to de¬ 
rive the Bogoliubov Hamiltonian) is spatially uniform. 
With these requirements, it can be shown that (r^Hk) 2 
is a Hermitian matrix where H is the BdG Hamiltonian 
introduced in Eq. (11) from the main text. This implies 
that the eigenvalues of {r^Hk ) 2 , namely (E^) 2 , will be 


real numbers and so will either be pure real or pure 
imaginary. This means that an imaginary part in the en¬ 
ergy spectrum will arise only in the region of overlap of 
particle and hole bands with vanishing real part of the en¬ 
ergy. When the zero value of the real energy spectrum is 
situated in the gapped region, no dynamical instability is 
induced except in the crossing edge states. In this sense, 
these symmetries protect the bulk bands from becoming 
dynamically unstable. While the model we consider in 
the main text possesses these symmetries, it is possible 
that less stringent conditions exist. 
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